In various fields, there have been attempts to analyze the behavior of an object by reproducing the motion of the object on a computer. The reproduction of object motion by a computer may be widely used to analyze, for example, the motion of a golf ball hit by a club, the resistance of a building to a tsunami, the behavior of a rubber film attached to a container for containing a liquid, etc.
Examples of objects whose motion is analyzed include fluids and elastic bodies. Numerical methods such as the finite difference method, the finite element method, and the finite volume method, which obtain an approximate solution of a differential equation based on the numerical mesh, have been widely used to solve a continuum of fluid or elastic body. In recent years, numerical methods have been advanced so that they can be used in fields such as computer aided engineering (CAE) and to solve a problem of interaction between a fluid and a structure. However, handling of calculation methods using the mesh is complicated when they are applied to, for example, a problem involving an interface such as a free surface or the problem of the interaction between a fluid and a structure involving a moving boundary, and therefore programming for such calculation methods tends to be difficult.
There also exist particle methods that analyze an object as a set of particles. Well-known particle methods include a Moving Particle Semi-implicit (MPS) method and a Smoothed Particle Hydrodynamics (SPH) method (see, for example, S. Koshizuka, A. Nobe and Y. Oka, “Numerical Analysis of Breaking Waves Using The Moving Particle Semi-Implicit Method”, International Journal for Numerical Method in Fluids, 26: 751-769 (1998); and J. J. Monaghan, “Smoothed Particle Hydrodynamics”, Annu. Rev. Astron. Astrophys. 30: 543-74 (1992)). Because particle methods do not require a special procedure for handling a moving boundary, they have become widely used even for the problem of the interaction between the fluid and the structure.
For example, there exists an analysis apparatus that records data regarding the velocity, position, and pressure for each of first particles representing a fluid and each of second particles representing a wall contacting the fluid, and analyze the behavior of the fluid taking into account a contact angle (see, for example, Japanese Laid-Open Patent Publication No. 2008-111675).
In calculations according to a particle method, a continuum is divided into particle distributions, and a resultant force of stress of the continuum is expressed by a sum of forces of interaction between particles. When the particles come close to each other, a repulsive force is generated as a result of increased pressure. For this reason, the particles do not become excessively dense or sparse.
Interaction is expressed using kernel functions that are used to construct a continuum field from particle distributions. In the SPH method, for example, a problem for solving the motion of a fluid is reduced to ordinary differential equations expressed by formulas (1) through (4) below. Formulas (1) and (2) discretize the conservation of momentum of fluid, and formula (3) discretizes the conservation of mass of fluid. Formula (4) is a simple example of a fluid equation where the density and the pressure have a linear relationship. Instead of formula (4), any other equation of state may be used. In the formulas, a subscript “i” indicates a particle number. The position vector, velocity vector, density, and pressure of an i-th particle are represented by #ri, #vi, ρi, and pi, respectively. Here, “#” indicates that the following alphabetic character represents a vector. Also, mj indicates the mass of a j-th particle.
                                          ⅆ                          r              i                                            ⅆ            t                          =                  v          i                                    (        1        )                                                      ⅆ                          v              i                                            ⅆ            t                          =                  -                                    ∑              j                                                                    ⁢                                                            m                  j                                ⁡                                  (                                                                                    p                        i                                            +                                              p                        j                                                                                                            ρ                        i                                            ⁢                                              ρ                        j                                                                              )                                            ⁢                                                ∂                                      W                    ⁡                                          (                                                                                                                                                            r                              i                                                        -                                                          r                              j                                                                                                                                ,                        h                                            )                                                                                        ∂                                      r                    i                                                                                                          (        2        )                                                      ⅆ                          ρ              i                                            ⅆ            t                          =                              ∑            j                                                          ⁢                                    m              j                        ⁢                                          ρ                i                                            ρ                j                                      ⁢                                          (                                                      v                    i                                    -                                      v                    j                                                  )                            ·                                                                    r                    i                                    -                                      r                    j                                                                                                                                r                      i                                        -                                          r                      j                                                                                                              ⁢                                          ∂                                  W                  ⁡                                      (                                                                                                                                                r                            i                                                    -                                                      r                            j                                                                                                                      ,                      h                                        )                                                                              ∂                                  (                                                                                                        r                        i                                            -                                              r                        j                                                                                                  )                                                                                        (        3        )                                          p          i                =                  c          ⁡                      (                                          ρ                i                            -                              ρ                0                                      )                                              (        4        )            
For “W” in formula (3), for example, a cubic spline function may be used as a kernel function. When the kernel function is expressed by formula (5) below, “h” indicates the radius of the support of the kernel function.
                              W          ⁡                      (                          r              ,              h                        )                          =                  {                                                                                          (                                          1                      -                                              1.5                        ⁢                                                                              (                                                          r                              h                                                        )                                                    2                                                                    +                                              0.75                        ⁢                                                                              (                                                          r                              h                                                        )                                                    3                                                                                      )                                    /                  β                                                                              0                  ≤                                      r                    h                                    <                  1                                                                                                      0.25                  ⁢                                                                                    (                                                  2                          -                                                      r                            h                                                                          )                                            3                                        /                    β                                                                                                1                  ≤                                      r                    h                                    <                  2                                                                                    0                                                              2                  ≤                                      r                    h                                                                                                          (        5        )            
In other words, “h” represents the radius of influence of particles, and is often set at a value that is two to three times greater than the initial average distance between particles. Also in formula (5), β indicates a value that is adjusted so that the total space integral of the kernel function becomes 1. When the kernel function is two-dimensional, β is set at 0.7 πh2. When the kernel function is three-dimensional, β is set at πh3.
In the SPH method, formula (3) indicates that density increases as particles come closer to each other, formula (4) indicates that pressure increases as the density increases, and formula (2) indicates that a force is generated in a direction from a higher pressure area toward a lower pressure area and results in a repulsive force between particles. FIG. 1 is a drawing illustrating a repulsive force f generated between particles i and j coming closer to each other.
Thus, in the related-art particle method, it is assumed that stationary particles are arranged on a boundary and a repulsive force is generated between the stationary particles and particles coming closer to the boundary to express impermeability (which indicates that an object does not penetrate through a boundary of another object). In the case of a fluid-structure coupled problem, the interaction between a fluid and a structure is expressed based on a pressure applied to boundary particles of the structure. FIG. 2 is a drawing illustrating a particle D1 of a fluid receiving repulsive forces f from particles E1, E2, . . . of a wall.